/*
 * topological-hex, a program to compute combinatorial hexahedral meshes.
 *
 * Copyright (C) <2018> <Université catholique de Louvain (UCL), Belgique>
 *
 * List of the contributors to the development, description and complete
 * License: see LICENSE file.
 *
 * This program (topological-hex) is free software:
 * you can redistribute it and/or modify it under the terms
 * of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program (see COPYING file).  If not,
 * see <http://www.gnu.org/licenses/>.
 */

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stdbool.h>
#include <limits.h>
#include <string.h>

#include <pthread.h>

#include "errors.h"
#include "types.h"
#include "utils.h"
#include "vertex_set.h"
#include "mesh.h"
#include "solver.h"
#include "dfs.h"
#include "shellable_dfs.h"
#include "parallel_dfs.h"
#include "lower_bound.h"
#include "compatibility.h"

static void on_error(const char *message, void *program_name) {
  fprintf(stderr, "%s: %s\n", (const char*)program_name, message);
}

typedef struct options {
  const char *filename;
  const char *symmetries_filename;
  const char *precomputed_filename;
  const char *lower_bounds_filename;

  vertex_index max_v_int;
  hex_index max_hex;

  uint8_t flags;

  solver_callback on_solution;
} options;

static
void parse_options(int argc, char **argv, options *opts);

static
void print_mesh(const solver *solver, uint32_t num_hexes, void *data);

static
void print_text(const solver *solver, uint32_t num_hexes, void *data);

static
void print_sexp(const solver *solver, uint32_t num_hexes, void *data);

int main(int argc, char **argv) {
  options opts;
  parse_options(argc, argv, &opts);

  set_error_log_fn(on_error, argv[0]);

  if (opts.lower_bounds_filename) {
     if (load_lower_bounds(opts.lower_bounds_filename) != SUCCESS)
      fprintf(stderr, "%s: failed to load lower bounds from file `%s'.\n",
              argv[0], opts.lower_bounds_filename);
  }

  if ((opts.flags & SOLVER_SHELLABLE_ONLY) && opts.precomputed_filename) {
    if (load_known_solutions(opts.precomputed_filename) != SUCCESS)
      fprintf(stderr, "%s: failed to load precomputed solutions from file `%s'.\n",
              argv[0], opts.precomputed_filename);
  }

  int status = 0;

  mesh mesh;
  if (mesh_load(&mesh, opts.filename) != SUCCESS) {
    status = 1;
    fprintf(stderr, "%s: failed to load mesh from file: %s\n", argv[0], opts.filename);
    goto fail_mesh_load;
  }

  if (opts.max_hex < mesh.num_quads / 6) opts.max_hex = mesh.num_quads / 6;

  if (mesh_reserve(&mesh, opts.max_hex) != SUCCESS) {
    status = 1;
    fprintf(stderr, "%s: failed to reserve memory for hexahedra\n", argv[0]);
    goto fail_mesh_reserve;
  }

  uint32_t num_v = mesh.num_vertices + opts.max_v_int;
  if (num_v >= VERTEX_SET_MAX) num_v = VERTEX_SET_MAX - 1;

  solver solver;
  if (solver_init(&solver, &mesh, num_v, opts.flags) != SUCCESS) {
    status = 1;
    fprintf(stderr, "%s: failed to allocate memory for the solver\n", argv[0]);
    goto fail_solver_init;
  }

  if (opts.symmetries_filename) {
    vertex_index *symmetries;
    vertex_index num_boundary_vertices;
    size_t num_symmetries;

    if (read_symmetries(opts.symmetries_filename, &symmetries, &num_boundary_vertices,
                        &num_symmetries) == SUCCESS) {
      solver_set_symmetries(&solver, num_symmetries, num_boundary_vertices,
                            symmetries);
    }
    else {
      fprintf(stderr, "%s: failed to load symmetries from file: `%s'\n",
              argv[0], opts.symmetries_filename);
      solver_compute_symmetries(&solver);
    }
  }
  else
    solver_compute_symmetries(&solver);

  pthread_mutex_t lock;
  if (pthread_mutex_init(&lock, NULL) != 0) {
    status = 1;
    fprintf(stderr, "%s: failed to initialize lock\n", argv[0]);
    goto fail_init_lock;
  }

  if (solver_run_parallel(&solver, opts.on_solution, &lock) != SUCCESS) {
    status = 1;
    fprintf(stderr, "%s: failed to run the solver in parallel\n", argv[0]);
    goto fail_parallel_solve;
  }
  
fail_parallel_solve:
  pthread_mutex_destroy(&lock);
fail_init_lock:
  solver_release(&solver);
fail_solver_init:
fail_mesh_reserve:
  mesh_release(&mesh);
fail_mesh_load:
  release_known_solutions();
  release_lower_bounds();
  return status;
}

static
void print_usage(const char *program_name, int status);

static
void print_version(void);

static
void parse_options(int argc, char **argv, options *opts) {
  opts->symmetries_filename = NULL;
  opts->precomputed_filename = NULL;
  opts->lower_bounds_filename = NULL;

  opts->max_hex = 64;
  opts->max_v_int = VERTEX_SET_MAX;

  opts->flags = SOLVER_BIPARTITE_GRAPH;

  opts->on_solution = print_sexp;

  for (int i = 1; i < argc; ) {
    if (strcmp(argv[i], "-H") == 0 || strcmp(argv[i], "--h-max") == 0) {
      if (i + 1 == argc) print_usage(argv[0], 1);

      char *end = NULL;
      opts->max_hex = strtoul(argv[i + 1], &end, 10);
      if (!end || *end != '\0') {
        fprintf(stderr, "%s: error: invalid limit to the number of hexahedra `%s'\n",
                argv[0], argv[i+1]);
        print_usage(argv[0], 1);
      }

      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));
      argc -= 2;
    }
    else if (strcmp(argv[i], "-V") == 0 || strcmp(argv[i], "--v-max") == 0) {
      if (i + 1 == argc) print_usage(argv[0], 1);

      char *end = NULL;
      opts->max_v_int = strtoul(argv[i + 1], &end, 10);
      if (!end || *end != '\0') {
        fprintf(stderr, "%s: error: invalid limit to the number of vertices `%s'\n",
                argv[0], argv[i+1]);
        print_usage(argv[0], 1);
      }

      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));
      argc -= 2;
    }
    else if (strcmp(argv[i], "-S") == 0 || strcmp(argv[i], "--symmetries") == 0) {
      if (i + 1 == argc) print_usage(argv[0], 1);
      opts->symmetries_filename = argv[i + 1];
      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));
      argc -= 2;
    }
    else if (strcmp(argv[i], "-s") == 0 || strcmp(argv[i], "--shellable-only") == 0) {
      opts->flags |= SOLVER_SHELLABLE_ONLY;
      memmove(argv + i, argv + i + 1, (argc - i - 1) * sizeof(*argv));
      argc--;
    }
    else if (strcmp(argv[i], "--non-bipartite-graph") == 0) {
      opts->flags &= ~SOLVER_BIPARTITE_GRAPH;
      memmove(argv + i, argv + i + 1, (argc - i - 1) * sizeof(*argv));
      argc--;
    }
    else if (strcmp(argv[i], "-p") == 0 || strcmp(argv[i], "--precomputed") == 0) {
      if (i + 1 == argc) print_usage(argv[0], 1);
      opts->precomputed_filename = argv[i + 1];
      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));
      argc -= 2;
    }
    else if (strcmp(argv[i], "-L") == 0 || strcmp(argv[i], "--lower-bounds") == 0) {
      if (i + 1 == argc) print_usage(argv[0], 1);
      opts->lower_bounds_filename = argv[i + 1];
      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));
      argc -= 2;
    }
    else if (strcmp(argv[i], "-f") == 0 || strcmp(argv[i], "--format") == 0) {
      if (i + 1 == argc) print_usage(argv[0], 1);

      if (strcmp(argv[i + 1], "sexp") == 0)
        opts->on_solution = print_sexp;
      else if (strcmp(argv[i + 1], "mesh") == 0)
        opts->on_solution = print_mesh;
      else if (strcmp(argv[i + 1], "text") == 0)
        opts->on_solution = print_text;
      else {
        fprintf(stderr, "%s: error: invalid output format `%s'\n",
                argv[0], argv[i+1]);
        print_usage(argv[0], 1);
      }

      memmove(argv + i, argv + i + 2, (argc - i - 2) * sizeof(*argv));
      argc -= 2;
    }
    else if (strcmp(argv[i], "--help") == 0 || strcmp(argv[i], "-h") == 0) {
      print_usage(argv[0], 0);
    }
    else if (strcmp(argv[i], "--version") == 0 || strcmp(argv[i], "-v") == 0) {
      print_version();
    }
    else
      i++;
  }

  if (argc != 2)
    print_usage(argv[0], 1);

  opts->filename = argv[1];
}

static
void print_usage(const char *program_name, int status) {
  FILE *out = status == 0 ? stdout : stderr;

  fprintf(out, "Usage: %s [OPTIONS] INPUT_MESH\n", program_name);
  fprintf(out, "\n");
  fprintf(out, "Lists the hexahedral meshes bounded by the quad mesh stored "
          "in INPUT_MESH.\n");
  fprintf(out, "\n");
  fprintf(out, "  -H, --h-max N\n");
  fprintf(out, "     Specifies the maximum number of hexahedra.\n");
  fprintf(out, "  -V, --v-max N\n");
  fprintf(out, "     Specifies the maximum number of interior vertices.\n");
  fprintf(out, "  -S, --symmetries FILENAME\n");
  fprintf(out, "     Explicitly specifies the symmetries that should be accounted for from\n"
               "     a file. Hexahedral meshes that are equivalent after applying one of\n"
               "     those symmetries are only output once. By default, all symmetries of\n"
               "     the input mesh are computed automatically.\n");
  fprintf(out, "  -s, --shelable-only\n");
  fprintf(out, "     Only list shellable hexahedral meshes. This can be significantly faster\n"
               "     if not all solutions are needed.\n");
  fprintf(out, "  --non-bipartite-graph\n");
  fprintf(out, "     Disable optimizations that assume the mesh graph to be bipartite. This is needed\n"
               "     to enumerate all solutions for certain non-ball domains (e.g. the torus).\n");
  fprintf(out, "  -p, --precomputed FILENAME\n");
  fprintf(out, "     Loads a list of hexahedral meshes from FILENAME, alonge with their boundaries.\n"
               "     A solution is output whenever the boundary of the unmeshed region matches\n"
               "     one of the meshes in this file. Duplicate solutions may appear.\n"
               "\n"
               "     Ignored unless the `--shelable-only' flag was also provided.\n");
  fprintf(out, "  -L, --lower-bounds FILENAME\n");
  fprintf(out, "     Loads a list of quad meshes along with the minimum number of hexahedra and \n"
               "     vertices needed to mesh them. Used to speed up the search by pruning nodes\n"
               "     when these quad meshes are encountered.\n");
  fprintf(out, "  -f, --format FORMAT\n");
  fprintf(out, "    sexp: One S-expression per mesh, separated by new lines.\n");
  fprintf(out, "    text: Human-readable format.\n");
  fprintf(out, "    mesh: MEDIT-compatible representation of solutions.\n");
  fprintf(out, "  -v, --version\n");
  fprintf(out, "     Prints information about the version of the software.\n");
  fprintf(out, "  -h, --help\n");
  fprintf(out, "     Displays this message.\n");
  fprintf(out, "\n");
  exit(status);
}

static
void print_version(void) {
  printf("topological-hex 0.3.0\n"
         "Copyright (C) 2019 Université catholique de Louvain (UCLouvain), Belgique\n"
         "License GPLv3+: GNU GPL version 3 or later <https://gnu.org/licenses/gpl.html>.\n"
         "This is free software: you are free to change and redistribute it.\n"
         "There is NO WARRANTY, to the extent permitted by law.\n"
         "\n"
         "Written by Kilian Verhetsel <kilian.verhetsel@uclouvain.be>.\n");
  exit(0);
}

static
void print_mesh(const solver *solver, uint32_t num_hexes, void *data) {
  pthread_mutex_t *lock = data;

  pthread_mutex_lock(lock);

  printf("Hexahedra\n");
  printf("%u\n", num_hexes);
  for (hex_index i = 0; i < num_hexes; i++) {
    for (hex_vertex_index j = 0; j < 8; j++) {
      printf("%u ", solver->hexes[8*i+j] + 1);
    }
    printf("0\n");
  }

  pthread_mutex_unlock(lock);
}

static
void print_text(const solver *solver, uint32_t num_hexes, void *data) {
  pthread_mutex_t *lock = data;

  pthread_mutex_lock(lock);

  printf("SOLUTION!\n");
  for (hex_index i = 0; i < num_hexes; i++) {
    printf("  %u) {", i);
    for (hex_vertex_index j = 0; j < 8; j++) {
      printf("%u", solver->hexes[8*i+j]);
      if (j + 1 != 8) printf(", ");
    }
    printf("}\n");

    for (hex_index j = i+1; j < num_hexes; j++) {
      if (!hex_compatibility(solver->hexes + 8*i, solver->hexes + 8*j)) {
        printf("  ERROR: Hex #%u and #%u are incompatible!\n", i, j);
      }
    }
  }

  pthread_mutex_unlock(lock);
}

static
void print_sexp(const solver *solver, uint32_t num_hexes, void *data) {
  pthread_mutex_t *lock = data;

  pthread_mutex_lock(lock);

  printf("#(");
  for (hex_index i = 0; i < num_hexes; i++) {
    printf("#(");
    for (hex_vertex_index j = 0; j < 8; j++) {
      printf("%u", solver->hexes[8*i+j]);
      if (j + 1 != 8) printf(" ");
    }
    printf(")");
    if ((uint32_t)(i + 1) != num_hexes) printf(" ");
  }
  printf(")\n");

  pthread_mutex_unlock(lock);
}
